%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% O-Ring Production Networks
%   Author: Demir, Fieler, Xu, Yang
%   Last updated: August 2022
%
% Purpose: Export Subsidy Counterfactual Analysis (Full Model)
%   - Resolve initial equilibrium
%   - Counterfactual 1: export subsidy (new baseline)
%   - Counterfactual 2: export subsidy w/ balanced trade
%   - Counterfactual 3: export subsidy w/ skill premium increase
%   - Counterfactual 4: export subsidy w/ agglomeration
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all; diary off; clc; clear; 

    
%% Resolve Initial Equilibrium
fprintf('\nInitial Equilibrium\n');

    % parameters
    define_param

    % add estiamted parameters 
    load('../../output/estimates/fullmodel_estimates.mat', 'est');     
    define_param_est     

    % data moment
    define_data_moment
    
    % constants
    define_omega_dist %firm type and density
    define_const %qgrid, phi_v, phy_y, D_F and other constants
    
    % solve initial equilibrium
    eqm = solve_eqm_outer(param, const);    
    define_wage_schedule %add wage schedule
    
    % calculate model moments and objfcn
    model = cal_model_moment(param, const, eqm); 
    
    % calculate vars for each firm type
    firm = cal_firm_vars(param, const, eqm);
   
  


%% Subsidy

    % subsidy: share t of the cost of finding customers in foreign
    const.t = 0.09;
    
    
%% Counterfactual #1: Export Subsidy As The New Baseline 
fprintf('\nExport Subsidy Counterfactual as the New Baseline\n');

    % solve the equilibrium
    eqm_cf1 = solve_eqm_cf1_outer(param, const, eqm);
    fprintf('n(q):\n'); fprintf('%.4f ', eqm_cf1.nq); fprintf('\n'); 

    % calculate vars for each firm type    
    firm_cf1 = cal_firm_vars(param, const, eqm_cf1);
    
    % counterfactual changes
    ge_cf1 = cal_chg_ge(eqm, firm, eqm_cf1, firm_cf1);       
   
    

%% Counterfactual #2: Export Subsidy w/ Balanced Trade    
fprintf('\nAlternative Counterfactual: Balanced Trade\n');

    % solve the equilibrium
    eqm_cf2 = solve_eqm_cf2_outer(param, const, eqm);
    fprintf('n(q):\n'); fprintf('%.4f ', eqm_cf2.nq); fprintf('\n'); 
   
    % calculate vars for each firm type    
    firm_cf2 = cal_firm_vars(param, const, eqm_cf2);
    
    % counterfactual changes
    ge_cf2 = cal_chg_ge(eqm, firm, eqm_cf2, firm_cf2);  
    

    
%% Counterfactual #3: Export Subsidy w/ Skill Premium Increase 
fprintf('\nAlternative Counterfactual: Increase Skill Premium\n');
           
    % increase efficiency wages
    %param.winc = 0.7933; %percentage increase at qmax
    param.winc = 0.8372; %percentage increase at qmax
    param.w = linspace(1, 1+param.winc/100, param.Q); 
    param.w_cf3 = param.w; %save   
    
    % redefine constant
    define_const            
    
    % solve the equilibrium
    eqm_cf3 = solve_eqm_cf1_outer(param, const, eqm);
    fprintf('n(q):\n'); fprintf('%.4f ', eqm_cf3.nq); fprintf('\n'); 
     
    % calculate vars for each firm type    
    firm_cf3 = cal_firm_vars(param, const, eqm_cf3);
    
    % counterfactual changes
    ge_cf3 = cal_chg_ge(eqm, firm, eqm_cf3, firm_cf3); 
    fprintf('calibrate winc s.t. percentage change in average quality for all firms is zero: %.6f\n', ge_cf3.overall_qlevel_all)


  
%% Counterfactual #4: Export Subsidy w/ Agglomeration 
fprintf('\nAlternative Counterfactual: Agglomeration\n');
    
    % reset w after previous counterfactual
    param.w = 1;
    
    % calculate employment changes in top quintile
    cutoff = 36+1;
    employment_q = (eqm.X0(cutoff:end)+eqm.X1(cutoff:end))./const.wage_grid(cutoff:end);
    employment_q_cf1 = (eqm_cf1.X0(cutoff:end)+eqm_cf1.X1(cutoff:end))./const.wage_grid(cutoff:end);
    employment_changes_cf1_q = employment_q_cf1./employment_q-1;
    employment_changes_cf1 = sum(employment_q_cf1)/sum(employment_q)-1

    sum(eqm.nq(cutoff:end))
    sum(eqm_cf1.nq(cutoff:end))

    % increase productivity at top quintile
    param.zinc = (1-param.alpha_m-param.alpha_s)*(0.25*(0.229+1/1.6)+0.75*(0.697-1/1.6)) ...
                 *employment_changes_cf1;
    const.omega0 = const.omega0 + param.zinc*(eqm.Qindicator(:,5)==1); 

    % redefine constant
    define_const        
    
    % solve the equilibrium
    eqm_cf4 = solve_eqm_cf1_outer(param, const, eqm);
    fprintf('n(q):\n'); fprintf('%.4f ', eqm_cf4.nq); fprintf('\n'); 
     
    % calculate vars for each firm type    
    firm_cf4 = cal_firm_vars(param, const, eqm_cf4);
    
    % counterfactual changes
    ge_cf4 = cal_chg_ge(eqm, firm, eqm_cf4, firm_cf4);  


    

%% Output tables and figures
diary off

    % create folder
    mkdir('../../output/figures') 
    mkdir('../../output/tables') 
    mkdir('../../output/temp') 

    % output table
    out_exportsubsidycf_tables

    % output figure: with baseline export shock counterfactual
    eqm_cf0 = eqm_cf1;
    clearvars -except param const eqm eqm_cf0
    save('../../output/temp/exportsubsidy_counterfactual.mat')

    run('../3 baseline counterfactual/main_baseline_counterfactual.m')
    clearvars -except eqm_cf1

    load('../../output/temp/exportsubsidy_counterfactual.mat')
    delete('../../output/temp/exportsubsidy_counterfactual.mat')

    out_exportsubsidycf_figures
    

    
